home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / ada / gwuada_9.zip / FARITH.C < prev    next >
C/C++ Source or Header  |  1993-07-27  |  13KB  |  509 lines
  1. /*
  2.  * Copyright (C) 1985-1992  New York University
  3.  * 
  4.  * This file is part of the Ada/Ed-C system.  See the Ada/Ed README file for
  5.  * warranty (none) and distribution info and also the GNU General Public
  6.  * License for more details.
  7.  
  8.  */
  9. /*
  10.  * multi-precision multiplication, division and exponentiation
  11.  * for fixed point computations.
  12.  */
  13.  
  14. #include <stdlib.h>
  15. #include "config.h"
  16. #include "int.h"
  17. #include "ivars.h"
  18. #include "intbp.h"
  19. #include "miscp.h"
  20. #include "farithp.h"
  21.  
  22. #define int_copy(u,v)   for(i=0;i<=FIX_LENGTH;i++) v[i] = u[i]
  23. #define int_conv(v,w)   w[0]=1, w[1]=v
  24.  
  25. static void int_norm(int *);
  26. static void int_print(int *);
  27.  
  28. #ifdef DEBUG_FARITH
  29. void int_print();
  30. #endif
  31.  
  32. /*
  33.  * M U L T I P L E     P R E C I S I O N     I N T E G E R
  34.  *
  35.  *       A R I T H M E T I C       P A C K A G E
  36.  *
  37.  *             Robert B. K. Dewar
  38.  *              June 16th, 1980
  39.  *
  40.  * This package of routines implements multiple precision integer
  41.  * arithmetic using what are called the "classical algorithms" in
  42.  * Knuth "Art of Programming", Volume 2, Section 4.3.1. The style
  43.  * of the algorithms follows Knuth fairly closely, and section
  44.  * 4.3.1 can be consulted for further theoretical details.
  45.  *
  46.  * Multiple precision integers are represented as tuples of small
  47.  * integers in the range 0 to FIX_BAS - 1, where FIX_BAS is a power of 2
  48.  * (actually 2 ** FIX_DIGS) which is the base used to represent the
  49.  * long integers. Essentially the successive elements of the tuple
  50.  * are the digits of the representation in base FIX_BAS. All integers
  51.  * are normalized so that the first digit is non-zero, except in
  52.  * the case of zero itself. All values are positive.
  53.  *
  54.  
  55.  * The constants FIX_BAS and FIX_DIGS must be defined as global constants
  56.  * in a program using these routines. It is assumed that the value
  57.  * FIX_BAS * FIX_BAS - 1 can be represented as an integer value.
  58.  
  59.  * The following routines are provided:
  60.  
  61.  *      int_norm        integer normalization
  62.  *    int_div        integer division
  63.  *    int_mul        integer multiplication
  64.  *    int_tom        conversion to multi-precision
  65.  *    int_tol        conversion to long integer
  66.  */
  67.  
  68.  
  69. void int_tom(int *v, long n)                                    /*;int_tom*/
  70. {
  71.     /* Convert a positive integer to a multiple precision integer. */
  72.  
  73.     int        d;
  74.  
  75.     /* Build up v in groups of FIX_BAS digits until n is depleted. */
  76.  
  77.     d = v[0] = 5;
  78.     while(d) {
  79.         v[d--] = n % FIX_BAS;
  80.         n /= FIX_BAS;
  81.     }
  82.     int_norm(v);
  83. }
  84.  
  85. void int_mp2(int *u, int p)    /*;int_mp2*/
  86. {
  87.     /* Multiplication by a power of 2. (Shift) */
  88.  
  89.     int     s;
  90.     int     m;
  91.     int     i,f,t,c;
  92.  
  93.     m = p % FIX_DIGS;
  94.     s = p / FIX_DIGS;
  95.     c = 0;
  96.  
  97.     if (m) { /* needs multiplication */
  98.         f = 2;
  99.         for (i=1;i<m;i++) f *= 2; /* f = 2 ** m */
  100.         if (f * u[1] > FIX_BAS) { /* carry: needs extension */
  101.             for (i = u[0]; i > 0; i--) {
  102.                 t = u[i] * f + c;
  103.                 u[i+1] = t % FIX_BAS;
  104.                 c = t / FIX_BAS;
  105.             }
  106.             u[1] = c;
  107.             u[0] += 1;
  108.         }
  109.         else {
  110.             for (i = u[0]; i > 0; i--) {
  111.                 t = u[i] * f + c;
  112.                 u[i] = t % FIX_BAS;
  113.                 c = t / FIX_BAS;
  114.             }
  115.             /* carry is zero now */
  116.         }
  117.     }
  118.     /* now we just have to add some zeros at the end */
  119.     if (s) { /* needs shift */
  120.         for (i=1;i<=s;i++)
  121.             u[u[0]+i] = 0;
  122.         u[0] += s;
  123.     }
  124. }
  125.  
  126. void int_mul(int *u, int *v, int *w)                                /*;int_mul*/
  127. {
  128.     /* Multiply unsigned integers, cf Knuth 4.3.1 Algorithm M
  129.      *
  130.      * Multiplication is similar to, but not identical with, the
  131.      * usual pencil and paper algorithm. The main difference is
  132.      * that the sum is accumulated as we go along, rather than
  133.      * forming all the partial sums and adding them up at the end.
  134.      *
  135.      */
  136.  
  137.     int        i, j, k, t;
  138.  
  139.     /* Initialize result to all zeroes(actually it is only absolutely
  140.      * necessary to initialize the last #v digits of w to zero).
  141.      */
  142.  
  143. #ifdef DEBUG_FARITH
  144.     printf("int_mul\n"); 
  145.     int_print(u); 
  146.     int_print(v);
  147. #endif
  148.     w[0] = u[0] + v[0];
  149.     for (i = 1; i <= FIX_LENGTH; i++) w[i] = 0;
  150.  
  151.     /* Outer loop, through digits of multiplier */
  152.  
  153.     for (j = v[0]; j > 0; j--) {
  154.         /* The inner loop, through the digits of the multiplicand, is
  155.          * similar to the inner loop of the addition, except that the
  156.          * product is added in, and the carry, k, can exceed 1.
  157.          */
  158.  
  159.         k = 0;
  160.         for (i = u[0]; i > 0; i--) {
  161.             t = u[i] * v[j] + w[i + j] + k;
  162.             w[i + j] = t % FIX_BAS;
  163.             k = t / FIX_BAS;
  164.         }
  165.  
  166.         /* The final step in the inner loop is to store the final
  167.          * carry in the next position in the result.
  168.          */
  169.  
  170.         w[j] = k;
  171.     }
  172.  
  173.     /* The result must be normalized(there could be one leading zero),
  174.      * and then the result sign attached to the returned value.
  175.      */
  176.  
  177.     int_norm(w);
  178. #ifdef DEBUG_FARITH
  179.     int_print(w);
  180. #endif
  181. }
  182.  
  183. long int_tol(int *t)                                            /*;int_tol*/
  184. {
  185.     /* Convert a multiple precision integer to a C long integer, if possible.
  186.      * Otherwise, return 'OVERFLOW'.
  187.      *
  188.      * First check if it overflows.
  189.      */
  190.  
  191.     long    x;
  192.     int        i;
  193.  
  194.     arith_overflow = 0;        /* reset overflow flag */
  195.  
  196.     if (t[0] > 5) {
  197.         arith_overflow = 1;
  198.         return 0;
  199.     }
  200.     if (t[0] == 5 && t[1] >= 8) { /* t >= 8 * 2**28 = 2**31 */
  201.         arith_overflow = 1;
  202.         return 0;
  203.     }
  204.  
  205.     x = t[1]; /* set first part of result */
  206.     for (i = 2; i <= t[0]; i++)
  207.         x = FIX_BAS * x + t[i];
  208.  
  209.     return (x);
  210. }
  211.  
  212. void int_div(int *u, int *v, int *q)                            /*;int_div*/
  213. {
  214.     /* Obtain quotient and remainder of signed integers,
  215.      * cf Knuth 4.3.1 Algorithm D
  216.      * Result is returned as a tuple [quotient,remainder].
  217.      *
  218.      * This is by far the most difficult of the four basic operations.
  219.      * This is because the paper and pencil algorithm involves certain
  220.      * amounts of guess work which cannot be programmed directly. The
  221.      * approach(which is analyzed in detail in Knuth) is to reduce
  222.      * the guess work by computing a rather good guess at each digit
  223.      * of the result, and then correcting if the guess is wrong.
  224.      *
  225.      * If the divisor is zero, return om.
  226.      */
  227.  
  228.     int        i, j, k, du, p, c, d;
  229.     int        rr, n, m, qe;
  230.  
  231.     /* A special case, if u is shorter than v, the result is 0 */
  232.  
  233.     if (u[0] < v[0]) {
  234.         int_conv(0,q);
  235.         return;
  236.     }
  237.  
  238.     /* The case of a one digit divisor is treated specially. Not only
  239.      * is this more efficient, but the general algorithm assumes that
  240.      * the divisor contains at least two digits.
  241.      */
  242.  
  243.     if (v[0] == 1) {
  244.         q[0] = u[0];
  245.  
  246.         /* Basically this case is straight-forward. Since we can represent
  247.          * numbers up to FIX_BAS * FIX_BAS - 1, we can do the steps of the
  248.          * division exactly without any need for guess work. The division is
  249.          *  then * done left to right (essentially it is the short division
  250.          * case).
  251.          */
  252.         rr = 0;
  253.         for (j = 1; j <= u[0]; j++) {
  254.             du = rr * FIX_BAS + u[j];
  255.             q[j] = du / v[1];
  256.             rr = du % v[1];
  257.         }
  258.     }
  259.     /* Here is where we must do the full long division algorithm */
  260.  
  261.     else {
  262.         n = v[0];
  263.         m = u[0] - v[0];
  264.         u[0] += 1; /* extend u */
  265.         for(i=u[0];i>1;i--) u[i] = u[i-1];
  266.         u[1] = 0;
  267.         q[0] = m+1;
  268.  
  269.         /* The first step is to multiply both the divisor and dividend
  270.          * by a scale factor. Obviously such scaling does not affect
  271.          * the quotient. The purpose of this scaling is to ensure that
  272.          * the first digit of the divisor is at least FIX_BAS / 2. This
  273.          * condition is required for proper operation of the quotient
  274.          * estimation algorithm used in the division loop. Note that
  275.          * we added an extra digit at the front of the dividend above.
  276.          */
  277.  
  278.         d = FIX_BAS /(v[1] + 1);
  279.  
  280.         c = 0;
  281.         for (j = u[0]; j > 0; j--) {
  282.             p = u[j] * d + c;
  283.             u[j] = p % FIX_BAS;
  284.             c = p / FIX_BAS;
  285.         }
  286.  
  287.         c = 0;
  288.         for (j = v[0]; j > 0; j--) {
  289.             p = v[j] * d + c;
  290.             v[j] = p % FIX_BAS;
  291.             c = p / FIX_BAS;
  292.         }
  293.  
  294.         /* This is the major loop, corresponding to long division steps */
  295.  
  296.         for (j = 1; j <=(m + 1); j++) {
  297.             /* Guess the next quotient digit by doing a division based on the
  298.              * leading digits. This estimate is never low and at most 2 high.
  299.              */
  300.             qe =(u[j] != v[1])
  301.               ? ((u[j] * FIX_BAS + u[j + 1]) / v[1]) :(FIX_BAS - 1);
  302.  
  303.             /* The following loop refines this guess so that it is almost always
  304.              * correct and is at worst one too high(see Knuth for proofs).
  305.              */
  306.  
  307.             while((v[2] * qe) >
  308.               ((u[j] * FIX_BAS + u[j + 1] - qe * v[1]) * FIX_BAS + u[j + 2])) {
  309.                 qe -= 1;
  310.             }
  311.  
  312.             /* Now(for the moment accepting the estimate as correct), we
  313.              * subtract the appropriate multiple of the divisor. This is
  314.              * similar to the inner loop of the mul